Equilibria and Stability Analysis

This chapter extends on our mathematical analysis of models by introducing the concepts of points of equilibria and stability analysis. These types of analyses allow you to determine many behaviors of a system without needing to fully solving its differential equation model.

Although the trajectory for the state variables in differential equation models generally cannot be determined analytically, several key properties of models can often still be determined. These properties include:

• The location of equilibrium points
• The stability of the equilibrium points

An equilibrium point is defined as a set of state variable values that will cause the system to cease changing. Once the system enters an equilibrium configuration, it will not leave that configuration without an external stimulus. For instance, in our exponential growth model a single equilibrium point exists: that of zero people. If the population is empty, then the population will not grow and instead remain at 0 indefinitely.

In the exponential growth population model there is only one equilibrium point (). In other models you may have multiple equilibrium points. In a model of a highly infectious, incurable disease you can imagine a system where two equilibrium points exist: one where no one is infected and a second point where everyone is infected. As long as there were no infectious individuals, the population would remain healthy. If just a single infected individual were introduced into the population, the infection would, however, spread until everyone was infected and the population would then remain at that point (remember this hypothetical disease is incurable).

Multiple types of equilibria exist. Figure 1 illustrates what is known as the stability of equilibrium points. Each of the three panes in this figure shows a different form of equilibrium for the ball. In all three the balls are in equilibrium; if the no external forces come into play, the balls will not move. What differs in each of the three is what occurs if the balls are displaced by a small amount.

Figure 1. Three different types of stability.

Stable Equilibrium : In this type of equilibrium the ball will return to its original position if it is displaced. The structure of the system is such that the system is naturally attracted to the point of equilibrium. To use the physical metaphor, the equilibrium is at the bottom of a dip and the system naturally rolls into it.

Unstable Equilibrium : Here the ball will move further and further away from the point of equilibrium if it is displaced by even a small amount. The equilibrium is unstable in that if we are just a small distance away from it, we move further away from it. To use the physical metaphor, the equilibrium is at the top of the hill and the system will move away from it unless it is placed at the exact point of equilibrium.

“Neutrally Stable Equilibrium” : This is a less common form of equilibrium and goes by several different names. In this case if the ball is moved it will stay fixed at its new location. It will not move closer to or further from the original equilibrium. Of the three types of equilibrium, this one is of less interest or relevance in practice.

In the case of the highly infectious disease model, an equilibrium of everyone being healthy would be classified as an unstable equilibrium. The equilibrium would persist as long as no one brought the disease into the population (someone would not just spontaneously become ill), but if as little as a single sick person entered the population, the population would move further and further away from the equilibrium point of everyone being healthy and would never naturally return to it.

On the other hand, the equilibrium point of everyone being sick is a stable equilibrium, as no one recovers from the disease on their own. Even if you introduced healthy people into a population of sick individuals – moving the population away from the equilibrium – they too would eventually become sick, restoring the population to the equilibrium of everyone being sick.

Equilibrium Points

Often, we can determine the equilibrium points for a system without fully needing to solve the trajectory for the state variables. Let’s implement the simple disease model we’ve been discussing. We’ll do so for both a differential equation model and a System Dynamics model, but we’ll rely on the differential equation version to do our analytic analysis.

One way to express the differential version of the model is to define two state variables: the number of healthy people () and the number of sick people (). The rate of infection between sick and healthy people can be made a function of the number of people in each category. Clearly, if there are no sick people the infection rate is 0; but, just as clearly, if everyone is already sick then the infection rate will also be zero. One workable differential equation model to implement this behavior is shown below:

This model uses a single parameter () to control the infection rate. is a non-zero positive value; the smaller is, the slower the infection will progress and vice versa. This notation illustrates one of the clumsier aspects of implementing stock and flow models using differential equations. The flow values between two stocks have to be repeated twice; once for each of the two connected state variables’ derivatives.

Incurable Disease

This model illustrates stable and unstable equilibria using the scenario of an incurable disease in a population.

Analytically, finding the equilibria for differential equation models is by-and-large straightforward. We simply need to harness the definition of an equilibrium point: an equilibrium point is one where the state variables are constant and unchanging. Since the derivatives represent changes in the state variables, this statement is equivalent to saying the derivatives for the model are 0 at equilibrium points.

Based on this, in order to find the equilibrium points we simply need to set the derivatives in our model to 0 and solve the resulting equations. For the disease model we get:

The initial conditions will determine the equilibrium but they do not affect the existence of the equilibria. Furthermore, the two equations we have set to 0 are equivalent⁷⁰. We can simplify these equations to:

Simple inspection reveals that this equation is true if and only if either , , or . Thus we have mathematically shown that our equilibria are either when everyone is sick or everyone is healthy (or there is no infection whatsoever). As we said earlier, this is a trivial conclusion for this model. However, for more complex models this type of analysis can be very useful and will often reveal that equilibria are functions of the different parameter values in the model. They may enable you to explicitly determine how the equilibria changes as the model configuration changes.

Let’s try a more complex example. Remember the predator-prey model from earlier? We had the following set of equations to simulate the relationship between a moose and wolf population:

Let’s determine the equilibrium values for this model. As before, we start by setting the derivatives to 0:

Solving this set of equations is more difficult than for the disease model. However, a little bit of algebra reveals two solutions. The first is when and (there are no animals at all), and the second is when and . This illustrates the dependency of the equilibrium location on the values of the model parameters.

The Phase Plane

When looking at model results we have been focused on time series plots and we have mainly been interested in the trajectory of the variables and stocks over time. For the mathematical analysis of differential equations, however, the primary graphical tool is not this time series plot; instead it is what is known as a phase plane plot.

Phase planes are almost like scatterplots. They show one of the state variables plotted against another of the state variables. A scatterplot could be used to show the path for these two variables over the course of a simulation. In the predator-prey model the results of a scatterplot of the wolf and moose population will be an ellipsoid. The two populations will cycle continuously. A phase plane plot is similar to this, but rather than just showing one of these cycles for a given simulation run, the phase plane shows the trajectories for all combinations of moose and wolf population sizes.

Figure 2. Predator-prey phase plane plot. The trajectory for a single set of initial conditions is highlighted in red.

Figure 2 illustrates a phase plane plot for the predator-prey system. The trajectory for one set of parameter and state variable values is highlighted in red and, as expected, we see a continual oscillation. We can also see the trajectories for all the other combinations of state variables. We see that the system will always oscillate and the size of this oscillation depends on the initial conditions for the state variables. This illustration provides us with a good deal of information in a single graphic; the phase plane plot is a great way to summarize the behavior of a system with two state variables.

Let’s quickly explore the phase plane plots for a simpler system. Take a system consisting of two state variables, both of which grow (or decay) exponentially.⁷¹ These state variables will be assumed to be independent from each other, so the value of one does not affect the value of the other:

Clearly, there is an equilibrium point for this model at and . There are four general types of behavior around this equilibrium. 1) when and , 2) when and , 3) when and , and 4) when and . The phase planes for each of the four cases are shown in Figure 3.

Figure 3. Phase planes for a simple two state variable exponential growth model.

From these plots we can visually determine how the stability of the equilibrium point at changes as we change and . When and , we have a stable equilibrium. In all other cases we have an unstable equilibrium.

Stability Analysis

Now that we have learned how to analytically determine the location of equilibrium points, we may want to determine what type of stability occurs at these equilibria. As we stated earlier, for the incurable disease model it is trivial to conclude that the state of everyone being healthy is unstable, while the state of everyone being sick is stable. In more complex models, it may be harder to draw conclusions, or the stability of an equilibrium point may change as a function of the model’s parameter values. Fortunately, there is a general way to determine the precise stability nature of the equilibrium points analytically.

The procedure to do this is relatively straightforward, but the theory behind it can be difficult to understand. The first key principle that must be understood is that of “linearization”. To get a feel for linearization, let’s take the curve in Figure 4. Clearly this curve is not linear. It has lots of bends and does not look at all like a line.

Figure 4. As we zoom in on a function it becomes more and more linear.

If we zoom in on any one part of the curve, however, the section we are zoomed in on starts to straighten out. If we keep zooming in, we will eventually reach a point where the section we are zoomed in on is effectively linear: basically a straight line. This is true for whatever part of the curve we zoom in on⁷². The more bendy parts of the curve will just take more zooming to convert them to a line.

We can conceptually do the same for the equilibrium points in our phase planes. Even if the trajectories of the state variables in the phase planes are very curvy, if we zoom in enough on the equilibrium points, the trajectories at a point will eventually become effectively linear. The simple, two-state variable exponential growth model we illustrated with phase planes above is an example of a fully linear model. If we zoom in sufficiently on the equilibrium points for most models, the phase planes for the zoomed-in version of the model will eventually start to look like one of these linear cases.

Mathematically, we apply linearization to an arbitrary model by first calculating what is called the Jacobian matrix of the model. The Jacobian matrix is the matrix of partial derivatives of each derivative in the model with respect to each of the state variables:

The Jacobian is a linear approximation of our (potentially) non-linear model derivatives. Let’s take the Jacobian matrix for the simple exponential growth model:

This is complicated so don’t worry if you don’t completely understand it! Once you have the Jacobian, you calculate what are known as the eigenvalues of the Jacobian at the equilibrium points. This is also a bit complicated, so if your head is starting to spin, just skip forward in this chapter!

Nonetheless, eigenvalues and their sibling eigenvectors are an interesting subject. Given a square matrix (a matrix where the number of rows equals the number of columns), an eigenvector is a vector which, when multiplied by the matrix, results in the original vector multiplied by some factor. This factor is known as an eigenvalue and is usually denoted as . Given a matrix , an eigenvalue with associated eigenvector ; the following equation will be true:

Let’s look at an example for a matrix:

What eigenvector and eigenvalue combinations satisfy this equation? It turns out there are two key ones:

Naturally, any multiple of an eigenvector will also be an eigenvector. For instance, in the case above, and are also eigenvectors of the matrix.

We can interpret eigenvectors geometrically. Looking at the matrix case, we can think of a vector as representing a coordinate in a two-dimensional plane: . When we multiply our matrix by the point, we transform the point into another point in the two-dimensional plane. Due to the properties of eigenvectors, we know that when we transform an eigenvector, the transformed point will just be a multiple of the original point. Thus, when a point on an eigenvector of a matrix is transformed by that matrix, it will move inward or outward from the origin along the line defined by the eigenvector.

We can now relate the concept of eigenvalues and eigenvectors to our differential equation models. Take a look back at the phase planes for the exponential model example. For each phase plane, there are at least two straight lines of trajectories. The x- axis and the y-axis are the locations of these trajectories. A system on the x- or y-axis will remain on that axis as it changes. This indicates that for this model, the eigenvectors are the two axes, as a system on either of them does not change direction as it develops. That’s the definition of an eigenvector.

For our purposes though, we do not really care about the actual direction or angle for these eigenvectors. Rather, we care about whether the state variables move inward or outward along these vectors. We can determine this from the eigenvalues of the Jacobian matrix. If the eigenvalue for an eigenvector is negative, the values move inward along that eigenvector; if the eigenvalue is positive, the values move outward.

These eigenvalues tell us all we need to know about the stability of the system. Returning to our illustration of stability as a ball on a hill, we can think of eigenvalues as being the slopes of the hill around the equilibrium point. If the eigenvalues are negative, the ground slopes down towards the equilibrium point, forming a cup (leading to a stable equilibrium). If the eigenvalues are positive, the ground slopes away from the equilibrium point, creating a hill (leading to an unstable equilibrium).

Eigenvalues can be calculated straightforwardly for a given Jacobian matrix. Briefly, for the Jacobian matrix , the eigenvalues are the values that satisfy the following equation, where is the matrix determinant and is the identity matrix.

We can do a quick example of calculating the eigenvalues for the Jacobian matrix we derived for our two-state variable exponential growth model.

That is a fair amount of work to do. It’s even more complicated if you have more than two state variables. However, once you have gone through the calculations and determined the linearized eigenvalues for your equilibrium points, you know everything you might want to know about the stability of the system.

In the exponential growth model we can see that when the eigenvalues are both negative we have a stable equilibrium (refer to the graphs we developed earlier), while if either one is positive (or they both are) we have an unstable equilibrium. This is logical, since if either one is positive it pushes the system away from the equilibrium, making it unstable. If they are both negative, they both push the system toward the equilibrium point. Visualize the ball sitting in the cup or on the hill.

Looking at it this way, we realize that all we need in order to understand the stability of an equilibrium point are the eigenvalues of the Jacobian at the equilibrium point. This is an incredibly powerful tool. It reduces the complex concept of stability into an analytical procedure that can be applied straightforwardly.

Let’s now look at some more examples.

First let’s take our simple disease model from earlier. If you recall, that model was:

First let’s calculate the Jacobian for this model. We take the partial derivatives of the two derivatives with respect to each of the two state variables to create a two-by-two matrix:

Next, we evaluate this Jacobian at one of our equilibrium points. Let’s choose the one where the (no one is sick) and (where is the population size) so everyone is healthy:

We can now find the eigenvalues for this matrix. Once we go through the math we get two eigenvalues: 0 and . What do these mean? Well, since one of the eigenvalues is positive, this indicates we have movement away from the equilibrium point along at least one of the eigenvectors. The other vector has no movement (0 as the eigenvalue), but this one positive value will ensure we have an unstable equilibrium. Again, think of the ball. The positive eigenvalue indicates the ground slopes downward from the equilibrium point so a ball balanced on top of this hill will be very unstable.

Now let’s do the second equilibrium - the one where and (everyone is sick). Let’s evaluate the Jacobian at this equilibrium:

Now let’s find the eigenvalues for this matrix. Once we go through the math we get two eigenvalues: this time 0 and . Again, the 0 eigenvalue can be ignored, as it does not cause growth or change. However, the second eigenvalue is negative, indicating the system moves toward the equilibrium point again. Look back at our exponential growth phase planes. Negative coefficients indicate trajectories towards the equilibrium (create a cup for the ball). Thus, this second equilibrium is a stable one.

It’s time to look at a more complex example; we’ll consider our predator-prey model. First we calculate the Jacobian matrix for this model:

Now that we have the Jacobian, we’ll evaluate it at the trivial equilibrium of and . The resulting matrix is:

The eigenvalues of this matrix are and . One of the eigenvectors approaches the equilibrium and the other moves away from it. This means we have an unstable equilibrium. This is actually good news, as it indicates that the two animal populations will not spontaneously go extinct.

Let’s now evaluate the more complex equilibrium point we identified earlier of and . First we calculate the Jacobian at this point:

When we calculate the eigenvalues for this point we obtain and . Here the indicates the imaginary number . That’s a little strange, so how do we interpret this? Imaginary numbers in the eigenvalues indicate oscillations in the phase planes, thus this result means we have oscillations around the point of equilibrium. Since we have no real component in the eigenvalues, there is neither attraction towards the point of equilibrium or repulsion away from it, so we have a stable oscillation around the equilibrium.

Of course we already knew that from our simulations, but this stability analysis allows us to mathematically determine this relationship. You can see this is a very powerful tool. The following table summarizes the types of eigenvalues that can be found for a system with two state variables and their associated stabilities. In this table, “damped” oscillations refers to a system that oscillates around a point of stability. Over time, the oscillations will “dampen” growing smaller and smaller in size until the system arrives at the point of stability.

A “saddle” point is a point where one eigenvalue is positive and the other one is negative. In this case, one eigenvalue pushes the system towards stability, and the other eigenvalue pushes the system away from stability. The net effect of this process is actually instability. Only a single eigenvalue pushing the system away from stability is enough to make the system unstable.

Analytical vs. Numerical Analysis

The majority of this ILE has been focused on the numerical analysis of models and the qualitative conclusions that can be drawn from these results. In this chapter we introduced analytical tools that can be used – for the most part – to analyze the same models we have presented elsewhere. Take a moment to reflect on these different forms of analysis and what each one can offer.

The great benefit of the analytical techniques is that they can provide precise answers to the general behavior of the system. Most of these same answers can also be determined numerically (e.g., running the simulation many times and exploring the results), but those answers will be less precise and definite. If you manually attempt to explore the parameter space of your model, it is possible that you could miss some set of parameter values that will give you unexpected behavior. An analytical analysis may be fully comprehensive and can guarantee the completeness of your conclusions.

A weakness of analytical methods is that your model must be solvable analytically. This means that you will probably need to keep your model from growing too complex in order to keep it analytically tractable. Also, some common functions such as IfThenElse logic can make analytical work much more difficult. Further, some models may simply be impossible to analyze analytically, and may in fact be very simple in practice. For example, any model containing both and in the same equation will be intractable to many forms of analysis.

We think both analytical and numerical work are very applicable in practice. We do worry, though, about some of the analytical models and work we see presented or published. Sometimes these models seem to us to be much too simple to adequately represent the system they are supposed to be modeling. True, analytically the results of the models appear elegant and clear. But if the model is too simple to be relevant, these results have little use and may actually be very misleading. We sometimes worry that a focus on analytical work leads modelers to prioritize analytical tractability over model utility.⁷³ We believe a focus on analytical results can lead to reductionist models with reduced practical utility, and we caution modelers against becoming too focused on elegant solutions at the expense of relevance. Where available, more realistic models are preferable, even if they require numerical solutions rather than overly simplistic analytically solvable ones.

Summary

In this chapter we have introduced equilibrium and stability analyses. These are powerful techniques that can be used to draw definitive conclusions about the behavior of a system. These conclusions can supplement your simulation work to generate a comprehensive analysis of a dynamic system.